********************************************************************************
* Modelling Nash Equilibrium for industrial gas market                         *
* Lei Zhang                                                                    *
* 09/26/2013                                                                   *
********************************************************************************

SETS
    i   our_facilities          /F1*F2/
    ii  competitor's_facilities /FC1*FC2/
    j   markets                 /M1*M4/
    k   products                /LOX, LNI, GOX, GNI/
    t   time_periods            /T1*T4/;

*------------------------------------------------------------------------------*

SCALAR
    CU      upper_bound_for_capacity_expansion /100.0/
    M       Big-M_method                       /100.0/;

PARAMETERS
    fa(k)   price_elasticity_a  /LOX -0.01, LNI -0.01, GOX -0.01, GNI -0.01/
    fb(k)   price_elasticity_b  /LOX 2.5, LNI 2.5, GOX 2.5, GNI 2.5/
    fac(k)  price_elasticity_ac /LOX -0.005, LNI -0.005, GOX -0.005, GNI -0.005/
    fbc(k)  price_elasticity_bc /LOX 3.1, LNI 3.1, GOX 3.1, GNI 3.1/;

TABLE alphaa(i,k,t) Fixed_cost_to_expand_capacity
        T1      T2      T3      T4
F1.LOX  10.0    9.0     8.0     7.0
F1.LNI  10.0    9.0     8.0     7.0
F1.GOX  11.0    10.0    9.0     8.0
F1.GNI  11.0    10.0    9.0     8.0
F2.LOX  10.0    9.0     8.0     7.0
F2.LNI  10.0    9.0     8.0     7.0
F2.GOX  11.0    10.0    9.0     8.0
F2.GNI  11.0    10.0    9.0     8.0;

TABLE betaa(i,k,t) Unit_investment_cost_to_expand_capacity
        T1      T2      T3      T4
F1.LOX  1.1     0.9     0.7     0.5
F1.LNI  1.1     0.9     0.7     0.5
F1.GOX  1.2     1.0     0.8     0.6
F1.GNI  1.2     1.0     0.8     0.6
F2.LOX  1.0     0.9     0.8     0.7
F2.LNI  1.0     0.9     0.8     0.7
F2.GOX  1.1     1.0     0.9     0.8
F2.GNI  1.1     1.0     0.9     0.8;

TABLE gammaa(i,k,t) Unit_production_cost
        T1      T2      T3      T4
F1.LOX  0.6     0.6     0.6     0.6
F1.LNI  0.6     0.6     0.6     0.6
F1.GOX  0.5     0.5     0.5     0.5
F1.GNI  0.5     0.5     0.5     0.5
F2.LOX  0.5     0.5     0.5     0.5
F2.LNI  0.5     0.5     0.5     0.5
F2.GOX  0.6     0.6     0.6     0.6
F2.GNI  0.6     0.6     0.6     0.6;

TABLE Tr(i,j,k) Unit_transportation_cost
        LOX     LNI     GOX     GNI
F1.M1   0.10    0.10    0.10    0.10
F1.M2   0.27    0.27    0.19    0.19
F1.M3   0.22    0.22    0.16    0.16
F1.M4   0.10    0.10    0.30    0.30
F2.M1   0.70    0.70    0.26    0.26
F2.M2   0.34    0.34    0.29    0.29
F2.M3   0.08    0.08    0.44    0.44
F2.M4   0.50    0.50    0.20    0.20;

TABLE D(j,k,t) Demand_of_market
        T1      T2      T3      T4
M1.LOX  15.0    16.5    18.2    20.0
M1.LNI  35.0    38.5    42.4    46.6
M1.GOX  15.0    16.5    18.2    20.0
M1.GNI  35.0    38.5    42.4    46.6
M2.LOX  10.0    11.0    12.1    13.3
M2.LNI  30.0    33.0    37.4    41.1
M2.GOX  10.0    11.0    12.1    13.3
M2.GNI  30.0    33.0    37.4    41.1
M3.LOX   5.0     5.5     6.1     6.7
M3.LNI  15.0    16.5    18.2    20.0
M3.GOX   5.0     5.5     6.1     6.7
M3.GNI  15.0    16.5    18.2    20.0
M4.LOX  10.0    11.0    12.1    13.3
M4.LNI  35.0    38.5    42.4    46.6
M4.GOX  10.0    11.0    12.1    13.3
M4.GNI  35.0    38.5    42.4    46.6;

TABLE CC(ii,k) Capacity_of_competitor's_facility
        LOX     LNI     GOX     GNI
FC1     30      50      20      40
FC2     30      50      20      40;

* Parameters:
* CU: upper bound for capacity expansion
* alphaa(i,k,t): Fixed_cost_to_expand_capacity
* betaa(i,k,t): Unit_investment_cost_to_expand_capacity
* gammaa(i,k,t): Unit_production_cost
* Tr(i,j,k): Unit_transportation_cost
* D(j,k,t): Demand_of_market
* CC(ii,k): Capacity_of_competitor's_facility

*------------------------------------------------------------------------------*

VARIABLES
    npv             Net_present_value
    
    lambdaa(j,k,t)  kkt1;
    
BINARY VARIABLES
    x(i,k,t)        Selection_of_capacity_expansion

    z1(i,k,t)       Big-M_for_KKT
    z2(ii,k,t)      Big-M_for_KKT;

POSITIVE VARIABLES
    y(i,j,k,t)      Amount_of_product_that_facility_sells_to_market
    yc(ii,j,k,t)    Amount_of_product_that_competitor's_facility_sells_to_market
    c(i,k,t)        Capacity_of_facility
    dc(i,k,t)       Capacity_expansion
    p(i,k,t)        price_of_products
    pc(ii,k,t)      competitor's_price_of_products
    
    muu1(i,k,t)     kkt2
    muu2(ii,k,t)    kkt2

    cost            Markets'_cost;

*------------------------------------------------------------------------------*

EQUATIONS
    obj             Objective_function_maximize_the_NPV
    dmd(i,k,t)      Demand_satisfaction_to_the_installed_capacity
    epd(i,k,t)      Investment_decision_in_capacity_expansion
    upb(i,k,t)      Upper_bound_of_capacity_expansion
    tdm(j,k,t)      Demand_satisfaction_for_all_markets
    dmdc(ii,k,t)    Demand_satisfaction_to_the_installed_competitor's_capacity
    pe(i,k,t)       Price_elasticity
    pec(ii,k,t)     Competitor's_price_elasticity

    kkt11(i,k,t)    kkt11
    kkt12(i,k,t)    kkt12
    kkt21(ii,k,t)   kkt21
    kkt22(ii,k,t)   kkt22
    kkt31(i,j,k,t)  kkt21
    kkt32(ii,j,k,t) kkt22

    mc              Markets'_cost;

obj..               npv =e= SUM((i,j,k,t), p(i,k,t)*y(i,j,k,t))
                        - SUM((i,k,t), alphaa(i,k,t)*x(i,k,t))
                        - SUM((i,k,t), betaa(i,k,t)*dc(i,k,t)) 
                        - SUM((i,k,t), gammaa(i,k,t)*c(i,k,t))  
                        - SUM((i,j,k,t), Tr(i,j,k)*y(i,j,k,t));

dmd(i,k,t)..        c(i,k,t) =g= SUM(j, y(i,j,k,t));
epd(i,k,t)..        c(i,k,t) =e= c(i,k,t-1)+dc(i,k,t);
upb(i,k,t)..        dc(i,k,t) =l= CU*x(i,k,t);
tdm(j,k,t)..        SUM(i, y(i,j,k,t))+SUM(ii, yc(ii,j,k,t)) =e= D(j,k,t);
dmdc(ii,k,t)..      CC(ii,k) =g= SUM(j, yc(ii,j,k,t));
pe(i,k,t)..         p(i,k,t) =e= fa(k)*c(i,k,t)+fb(k);
pec(ii,k,t)..       pc(ii,k,t) =e= fac(k)*CC(ii,k)+fbc(k);

kkt11(i,k,t)..      muu1(i,k,t) =l= M*(1-z1(i,k,t));
kkt12(i,k,t)..      SUM(j, y(i,j,k,t))-c(i,k,t) =g= -M*z1(i,k,t);
kkt21(ii,k,t)..     muu2(ii,k,t) =l= M*(1-z2(ii,k,t));
kkt22(ii,k,t)..     SUM(j, yc(ii,j,k,t))-CC(ii,k) =g= -M*z2(ii,k,t);
kkt31(i,j,k,t)..    p(i,k,t)+muu1(i,k,t)+lambdaa(j,k,t) =e= 0.0;
kkt32(ii,j,k,t)..   pc(ii,k,t)+muu2(ii,k,t)+lambdaa(j,k,t) =e= 0.0;

mc..                cost =e= SUM((i,j,k,t), p(i,k,t)*y(i,j,k,t))
                        +SUM((ii,j,k,t), pc(ii,k,t)*yc(ii,j,k,t));

*------------------------------------------------------------------------------*

MODEL NASH_COMPETITOR /ALL/;
SOLVE NASH_COMPETITOR USING MINLP MAXIMIZING npv;

********************************************************************************
